function [VAR, Y, eff] = q3( N, n )
%Q3 Summary of this function goes here
%  This function calculates the variance of the MC Control Variate Method
%  Returns
%  VAR: vector of numerical variance of estimators
%  Y:   vector of MC estimation results: Crude, CV1, CV2

if nargin < 1
    N = 10000; n = 100;
elseif nargin < 2
    n = 100;  % in sample for b estimation
end
 

% (a) g(x) = x
%     y = [e^x - 1] / [e - 1]

% generate N + n uniform u and control variate x
x = rand(n,1);
y = (exp(x) - 1)/(exp(1) - 1);  % Generate in sample Yi
miu_x = 0.5;
cov_xy = cov(x,y);
b = polyfit(x,y,1); % Get b_hat as best linear fit slope
b_hat = b(1);

% Now out of sample
x = rand(N, 1);
y = (exp(x) - 1)/(exp(1) - 1);
y_cv = y - b_hat.*(x - miu_x);

var0 = var(y);
var1 = var(y_cv);

% MC Estimation
y0 = mean(y);
y1 = mean(y_cv);

% Now g(x) = 0.51x + 0.49x^2
clear x y miu_x cov_xy b_hat y_cv;
% Generate in sample Yi
x = rand(n,1);
y = (exp(x) - 1)/(exp(1) - 1);
x = 0.51*x + 0.49*x.^2;
miu_x = 0.51*0.5 + 0.49/3;
% Get b_hat
cov_xy = cov(x,y);
b_hat = cov_xy(1,2)/var(x);

% Now out of sample
x = rand(N,1);
y = (exp(x) - 1)/(exp(1) - 1);
x = 0.51*x + 0.49*x.^2;
y_cv = y - b_hat.*(x - miu_x);
var2 = var(y_cv);

% MC Estimation
y2 = mean(y_cv);

% Return the MC estimation results for Crude and CV methods
VAR = [var0, var1, var2];
Y = [y0, y1, y2];
eff = [var0/var1, var0/var2];

end

